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Large temperature and pressure fluctuations have a profound effect on 
turbulence development in transonic and supersonic jets. For high-speed, high- 
temperature jet flows, standard turbulence models lack the ability to predict the 
observed mixing rate of a shear layer. Several proposals to address this deficiency 
have been advanced in the literature to modify the turbulence transport equations 
in a variety of ways. In the present study, some of the most proven and simple 
modifications to two-equation turbulence models have been selected and 
implemented in NASA’s USM3D tetrahedral Navier-Stokes flow solver. The 
modifications include the addition of compressibility correction and pressure 
dilatation terms in the turbulence transport equations for high-speed flows, and the 
addition of a simple modification to the Boussinesq’s closure model coefficient for 
high-temperature jets. The efficacy of the extended models is demonstrated by 
comparison with experimental data for two supersonic axisymmetric jet test cases at 
design pressure ratio. 
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Nomenclature 

= freestream speed of sound (dimensional) 

= distance to the closest wall normalized by L 
= nozzle exit diameter normalized by L 
= turbulence kinetic energy normalized by a £ 

= reference length (dimensional) 

= freestream Mach number 
= stagnation pressure normalized by p^a^ 

= nozzle plenum stagnation pressure normalized by p^a^ 
= static pressure normalized by p^a^ 

= Reynolds number per unit length 

= nozzle exit radius normalized by L 
= freestream temperature (dimensional) 

= stagnation temperature normalized by T ^ 

= nozzle plenum stagnation temperature normalized by T 
= time normalized by L/ 

= Cartesian velocity components normalized by a 
= Cartesian coordinates normalized by L 
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= turbulence dissipation normalized by a^/v ^ 

= frees tream laminar kinematic viscosity (dimensional) 

= laminar kinematic viscosity normalized by v / 

= eddy viscosity normalized by v z 

= primary unknown in Spalart- Allmaras turbulence model, normalized by v ^ 
= frees tream density (dimensional) 

= density normalized by 

= stagnation energy per unit volume normalized by p^a^ 

= turbulence specific dissipation rate normalized by a^/vj 


I. Introduction 

T HE NASA Constellation program is relying heavily on Computational Fluid Dynamic (CFD) tools to 
complement ground-based wind tunnel testing to aid the design and aerodynamic database 
development for the ARES I Crew Eaunch Vehicle (CLV), Orion Crew Exploration Vehicle (CEV) and 
Launch Abort Vehicle (LAV). These next-generation space vehicles are currently under development to 
replace the Space Shuttle in the next decade. 

The ARES-series of launch vehicles utilize large numbers of control jets during ascent and stage 
separations to maintain stable flight. These control jets operate in close proximity to aerodynamic surfaces 
at very high temperatures and over a wide range of Mach numbers, from subsonic to supersonic. Similarly, 
powerful Abort Motor jets and smaller Abort Control Motor jets will dominate the aerodynamic 
characteristics of the LAV during a launch abort, should such an unfortunate emergency occur. As might be 
expected, maintaining accurate simulation of jet-dominated aerodynamics on these vehicles is a formidable 
task. Considerable effort is being expended to develop suitable experimental techniques and numerical 
models for simulating these flows. The TetrUSS USM3D flow solver [1] serves as the lead CFD tool 
within the ARES program to develop the computational aerodynamic database and loads for the vehicle up 
through supersonic speeds. The USM3D solver is also lending critical support to the aerodynamic 
development of the CEV and LAV vehicles. 

This paper focuses on the extension of the advanced numerical turbulence models within the USM3D 
Navier-Stokes flow solver to improve the simulation accuracy for jet-dominated flows. Within such flows, 
large temperature and pressure fluctuations have a profound effect on the development of turbulence for 
transonic and supersonic jets. For high-speed, high-temperature jet flows, standard turbulence models lack 
the ability to predict the observed mixing rate of a shear layer. Several proposals to address this deficiency 
have been advanced in the literature to modify the turbulence transport equations in a variety of ways. 
Some of the most proven and simple modifications to two-equation turbulence models have been selected 
and implemented in the USM3D flow solver. The modifications include the addition of compressibility 
correction and pressure dilatation terms in the turbulence transport equations for high-speed flows, and the 
addition of a simple modification to the Boussinesq’s closure model coefficient for high- temperature jets. 
The efficacy of the extensions is demonstrated by comparison with experimental data for two supersonic 
axisymmetric jet test cases at design pressure ratio. The new turbulence modeling capability is currently in 
heavy use within the NASA Constellation program following an extensive study on the Orion LAV with 
Abort Motors to determine which turbulence models best correlate with subsonic/transonic wind-tunnel 
data. 


II. Description of Methodology 

The NASA Tetrahedral Unstructured Software System (TetrUSS) [1] is a complete flow analysis 
system that has been widely used in industry and government since 1991. Its strength is derived from 
harnessing various component technologies into a user-friendly system to provide rapid, higher-order 
analysis and design capability to the applied aerodynamicists. The system consists of loosely integrated, 
user-friendly software that comprises a geometry setup utility GridTool [2], a tetrahedral grid generator 
VGRID [3], [4], a flow solver USM3D [5], [6], and post-processing visualization and data extraction 
utilities. The extensions addressed in this paper pertain to the USM3D flow solver, and will be described 
below. 
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A. USM3D Flow Solver 

USM3D is a tetrahedral cell-centered, finite volume Euler and Navier- Stokes (N-S) flow solver. 
Inviscid flux quantities are computed across each cell face using Roe’s [7] flux- difference splitting (FDS), 
Liou’s [8] Advection Upstream Splitting Method (AUSM+) or Toro’s HLLC scheme [9]. Spatial 
discretization is accomplished by a novel reconstruction process [10], based on an analytical formulation 
for computing solution gradients within tetrahedral cells. The solution is advanced to a steady-state 
condition by an implicit backward-Euler time-stepping scheme [11]. 

USM3D has several closure models for capturing flow turbulence effects. First is the Spalart-Allmaras 
(SA) one-equation model [12], which can be coupled with a wall function boundary condition to reduce the 
number of cells in the sublayer region of the boundary layer. However, the wall function is typically used 
only when extensive curvature-based separation is not expected in a flow simulation. The second model is 
the two-equation k-e turbulence model [13]. The third model is the Menter SST two- equation model [14]. 
These turbulence models have undergone rigorous testing on propulsion afterbody flows in Ref. [15]. In 
addition, nonlinear Algebraic Reynolds Stress Models (ARSM) of Shih et al. [16] and Girimaji [17] have 
also been incorporated in USM3D. The two algebraic Reynolds stress models give inherently better results 
than the linear stress model because of the explicit modeling of effects such as relaxation, and the specific 
inclusion of nonlinear anisotropic effects from the mean flow strain and vorticity. With a nonlinear model, 
the calculation of six independent, realizable Reynolds stress terms is possible. This type of detail is 
important for simulating complicated multidimensional flows. A capability to trip the flow at specified 
locations on aerodynamic surfaces has been implemented for the k-e turbulence model [18]. 

USM3D supports an array of useful boundary conditions (BCs). It contains the standard BCs of flow 
tangency or no-slip on solid surfaces, characteristic inflow/outflow for subsonic boundaries, and freestream 
inflow and extrapolation outflow for supersonic flow. It also contains some additional special BCs for jet 
exhaust and intake, a propeller/rotor actuator disk model, and passive porosity [19]. The most recent 
enhancements to USM3D are second-order time accuracy, Detached Eddy Simulation (DES) capability, 
dynamic grid motion [20], and overset-grid moving body capability [21]. 

The USM3D runs on massively parallel computers and clusters of personal computers. A grid- 
partitioning file is quickly generated during the preprocessing step. A global restart file is generated and 
saved during and after each run. 


B. Turbulence Model Enhancements 

The two-equation k-e and SST turbulence models have been modified for jet flow computations to 
improve the prediction of jet spreading rates. Compressibility correction formulations of Sarkar [22], and 
Wilcox [23], and pressure dilatation terms have been added to both two-equation models following the 
approach of Suzen and Hoffman [24]. In addition, a simple modification to the Boussinesq’s closure model 
coefficient has been implemented in the k-e turbulence model for high- temperature jets [25]. 


7. Mean Flow Governing Equations 


Reynolds Averaged Navier- Stokes (RANS) equations are the governing equations solved in this study. 
The perfect gas law is chosen to represent air properties. The mass, momentum, and energy conservation 
equations in the time-averaged conservative simplified form can be written as follows: 
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In the above equations, dy is Kronecker delta, Xy is the laminar shear stress, and y is the ratio of specific 
heats. The term u t u j is the turbulent shear stress, which will be defined later. 
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Three turbulence models are used in the current study: the standard SA model [12], linear k-e model 
[13] enhanced with compressibility corrections, pressure dilatation term and temperature correction as well 
as Menter SST model [14] enhanced with compressibility corrections and pressure dilatation term. The 
eddy-viscosity concept is used to model the Reynolds stresses. The algebraic Reynolds stress models of 
Shih-Zhu-Lumley (SZL) [16] and Girimaji [17] have been found to be the most beneficial, primarily in 
wall-bounded flows. Past applications of these models for jet flows using PAB3D code [25] have yielded 
unsteady solutions, and therefore these models have not been used in the current studies based on the k-e 
model. The standard SA model is included for reference, since it is a widely used turbulence model for 
computational applications. The governing equations for the three turbulence models are presented below. 

2. Standard SA Model 
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£2 is the vorticity magnitude. 
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As seen from the equations above, the trip term is not included in the USM3D implementation of the 
SA model. The constants used in the above model are defined below. 

c bl =0.1355, C„ 2 = 0.622, a -2/3, k = 0.41, C w> =C bi / K 2 + (l + C bi )/a 

C W2 = 0.3, C W} = 2, C v ,= 7.1, C h =1.1, C, 4 =2 

v wall = 0 is applied for the sub-layer resolving grid. The SA model can also be optionally coupled with 
the wall function formulation inside USM3D [6]. 
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3. k-e Model 
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In the above equations, T is the compressibility correction term and PD is the pressure dilatation term. 
T = a x M 2 for Sarkar [22] model, and T = a x (M 2 - M*)ff (AT, - M, o ) for Wilcox [23] model. 
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H(x) is Heaviside step function. H(x ) = 0 for M t ^ M and //(x) = 3/2 for M t > M tQ , M t ^ = 0.25 
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cq = 1.0, a 2 = 0.4, a 3 = 0.2 following Ref. [24]. 

/ 2 = l.-0.3exp(-R, 2 ), P,=T 
The boundary conditions for k and £ at the wall are: 
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C* is 0.09 for the linear model and is a function of vorticity, as well as strain tensors, for the nonlinear 
models. For the temperature-corrected linear model, C* is a function of stagnation temperature gradient; 
the exact functional relationship will be described later. 


The turbulent stress components for the linear model are formulated as: 
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In the SZL nonlinear model [16], the turbulent stresses are given by: 
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Girimaji nonlinear model [17] is given by: 


„ „ 2 „ . „ * fc 3 / M ^ 


Re, 


2v,S & - [C 2 W lVt S*,. - S tt W#) + G,(S ik S kj -{s^S^d,,)} 


®'V Re L/ 


W I 

7 2 ax , ax, 


s, = - 

* 7 2 


dU: du • 

L+ :L 

ax ax, 


The variables such as /3,C*,G 2 ,G 3 have been defined in the Appendix. In the temperature- corrected 
model [25], C* is prescribed as, 
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In the above formula, T t is stagnation temperature. 


4. SSTModel 
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The constants a^,a w ,y, and fi are evaluated using the general formula 0 = + (l- Fj)0 2 . The 

parameters 0i and 02 represent the values in the inner- and outer-region, respectively. The inner- and outer- 
set constants have been defined below. 
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III. Results and Discussion 


A. Eggers Nozzle 

The initial assessment of the present method is performed on a fully expanded supersonic cold free jet 
case that was experimentally investigated by Eggers [26]. This study involved a Mach 2.22 axisymmetric 
nozzle exhausting to quiescent air at the design pressure ratio. The total temperature of the jet was the same 
as the ambient temperature of 525 °R. For the numerical investigation, a fully unstructured tetrahedral grid 
was generated with normal-direction grid spacing equivalent to an average cell-center y + = 0.20 at the 
nozzle wall. The grid was generated inside one quadrant of the nozzle, and the symmetry-plane boundary 
condition was applied on two planes defining the quadrant. The tetrahedral grid consists of 1,313,752 cells, 
242,567 nodes, and 82,106 boundary triangles. A view of the surface grid on one of the symmetry planes is 
presented in Fig. 1. 

The computed solutions for this case were obtained by considering a freestream Mach number of 0.05 
to avoid the solution instability and convergence difficulty encountered by a compressible flow solver 
USM3D at very low speed. All three turbulence models, namely, Spalart-Allmaras (SA) one-equation 
model, Faunder-Sharma k-e linear model, and Menter SST model were applied in various computed 
solutions. The SA model was applied in its standard form without incorporating any modifications to either 
the equation itself or to the coefficients used in the equation. Both of the two-equation models were applied 
in their standard forms as well as in their modified forms that include pressure dilatation term along with 
compressible dissipation formulation of either Sarkar or Wilcox. For the sake of brevity, the enhanced 
model using Sarkar formulation will be denoted as “Sarkar_PD”, and that using Wilcox formulation will be 
denoted as “Wilcox_PD” in the remainder of this section. Figure 2 presents the convergence of mean-flow 
residual errors for all of the computed solutions obtained for this case. The SA model converges 
monotonically whereas k-e as well as SST models exhibit stalled convergence after about 3.5 orders of 
magnitude reduction in residual errors. Symmetry plane variations of Mach number, as obtained from the 
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SA, k-e Sarkar_PD, and Menter SST Sarkar_PD models, are depicted in Fig. 3. It is evident from this figure 
that the k-e Sarkar_PD model predicts a slightly longer potential core region as compared to that from the 
SA and SST Sarkar_PD models. 

Figures 4 and 5 present a detailed comparison of the computed and measured velocity profiles along the 
nozzle centerline, as well as along the radial direction at several axial locations. The standard SA model 
result is included here for reference, since it is a widely used turbulence model for computational 
applications. In Fig. 4, predicted velocity profiles from the standard k-e , k-e Sarkar_PD, and k-e 
Wilcox_PD models are compared with the corresponding SA model predictions and the measurements of 
Eggers. Figure 5 presents a comparison of the predicted velocity profiles from the standard SST, SST 
Sarkar_PD, SST Wilcox_PD, and SA models, as well as the measurements of Eggers. It is evident from 
both figures that the addition of pressure dilatation term and compressible dissipation correction markedly 
improves the accuracy of the solutions based on k-e and Menter SST models relative to the standard 
versions. Amongst the two formulations for compressible dissipation, the k-e model yields slightly better 
accuracy with the Wilcox formulation in the vicinity of the centerline. However, Sarkar formulation is in 
closer agreement to the measurements far away from the centerline. The SST model performs better with 
the Sarkar formulation. It is also seen that the predictions of radial spread of jet are in quite good agreement 
with the measurements up to the axial locations of 75 times the nozzle exit radius. It should be recalled that 
due to the limitation imposed by the compressible flow solver, the computations have been performed 
assuming a freestream Mach of 0.05 as opposed to using quiescent air condition of the experiment. This 
discrepancy may have been responsible for the poor agreement of the computations and measurements at 
the axial station located atx/r exit = 121.30. 

B. Seiner Nozzle 

This axisymmetric supersonic jet test case was experimentally investigated by Seiner [27]. The nozzle 
was operated at the design pressure condition over a wide range of plenum temperature. The nozzle exit 
Mach number was 2, and plenum temperature ranged from 563 °R to 2459 °R. This test case has been 
selected to assess the temperature correction formulation presently implemented in the k-e turbulence 
model. Computations have been performed for the nozzle plenum temperatures of 563 °R, 1359 °R and 
2009 °R. Computed results for the plenum temperature of 563 °R will be shown first to assess the impact, if 
any, of temperature correction formulation on the cold jet. The results for the plenum temperature of 2009 
°R will be presented next, as this specific condition was used in the calibration of the temperature 
correction model in Ref. [25]. This study will establish that no further tuning of the model is warranted due 
to different solution methodology employed in various CFD codes. Finally, the computational results will 
be presented for an intermediate plenum temperature of 1359 °R to assess whether the model remains 
effective at a plenum temperature different from the one used for the model calibration. Only the k-e 
turbulence model has been used in the computations. However, several different variants of k-e models 
such as the standard formulation, Sarkar_PD, and Wilcox_PD formulations - with and without temperature 
correction term - have been applied. 

A fully unstructured tetrahedral grid was generated inside one quadrant of the nozzle and symmetry- 
plane boundary condition was applied on two planes defining the quadrant. The normal-direction grid 
spacing equivalent to an average cell-center y + = 0.12 was applied at the nozzle wall. The tetrahedral grid 
consists of 2,706,755 cells, 496,241 nodes, and 157,172 boundary triangles. Identical grid has been used for 
computing the cold as well as hot jet solutions. Although not shown, all solutions converged to a level of 
about 3.5 orders of reduction in the mean-flow residual errors. 

Figure 6 shows the comparison of the computed and measured profiles of Mach number and stagnation 
pressure along the nozzle centerline for the cold jet case. It is evident from the figure that the solutions 
obtained with and without temperature correction model are almost identical. This is only to be expected as 
the total temperature of nozzle plenum is identical to that of ambient air, and the flow should be devoid of 
total temperature gradients to activate the temperature correction model. It is observed from the figure that 
the Sarkar_PD formulation compares very well with the measured data. 

Next, the accuracy of the present temperature correction formulation is assessed via two hot jet cases. 
Initially, nozzle plenum stagnation temperature of 2009 °R has been used for this purpose. Figure 7a 
presents symmetry-plane variation of the normalized stagnation temperature. It can be seen that the 
stagnation temperature gradient locates itself in the shear layer region. In the present temperature correction 
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formulation, this gradient is used for increasing the C * in the shear layer region as shown in Fig. 7b. The 
increased C* results in the increased eddy viscosity that produces enhanced mixing in the shear layer. 

Figure 8 presents a comparison of the computed and measured centerline variations of Mach number, 
stagnation pressure, and stagnation temperature for the nozzle plenum temperature of 2009 °R. It is clear 
from this figure that the temperature correction formulation in conjunction with the Sarkar_PD model 
provides a very good comparison with the measured data. Further assessment of the temperature correction 
model is performed using an intermediate temperature of 1359 °R for the nozzle plenum. Figure 9 depicts 
favorable comparisons of the centerline variations of Mach number, stagnation pressure, and stagnation 
temperature obtained from the temperature- corrected Sarkar_PD computations and the experimental 
measurements. It can be concluded from Fig. 6, 8 and 9 that the temperature correction model provides 
reasonable and consistent correlation with the measured data across the temperature range, from ambient to 
2009 °R. 

IV. Concluding Remarks 

The two-equation turbulence models in the TetrUSS USM3D Navier-Stokes flow solver have been 
extended to account for compressibility and temperature effects. These extensions were implemented in 
response to a need in the NASA Constellation program to improve the aerodynamic and loads prediction 
capability of the solver for jet-dominated flows that are often encountered on the ARES Crew Launch 
Vehicle and Orion Launch Abort Vehicle. 

The standard forms of k-e and SST two-equation turbulence models are enhanced with selected 
published compressibility correction and pressure dilatation terms, and assessed on two supersonic cold jet 
cases at design pressure conditions. These studies have shown that while the standard form of two-equation 
model offers no apparent advantage over the SA one-equation model, incorporation of the compressibility 
correction and the pressure dilatation terms significantly improves the accuracy of jet spreading rates 
predicted by the two- equation models. 

A temperature correction formulation is also implemented in the k-e model to improve the 
computational accuracy for hot jets. Good correlation with experimental data has been demonstrated for 
three different temperatures of nozzle plenum for a supersonic jet case at the design pressure condition. 

Additional work is underway to evaluate the new turbulence model extensions for subsonic/transonic jet 
flows. Meanwhile, this new capability is currently in heavy use within the NASA Constellation program. 

Appendix 

The following functions and variables are used in the SZL [16] algebraic Reynolds stress model: 
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For the Girimaji [17] model: 
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Figure 1. Symmetry-plane view of unstructured tetrahedral grid for Eggers nozzle. 




(a) k-e model 


(b) SST model 


Figure 2. Convergence history of mean flow residuals for SA model and (a) k-e model (b) SST 
model. 
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(c) SST Sarkar PD model 


Figure 3. Contours of computed Mach number on Eggers nozzle symmetry plane for (a) SA model 
(b) k-e model with Sarkar compressibility correction and pressure dilatation (PD) term, (c) SST 
model with Sarkar compressibility correction and pressure dilatation (PD) term. 
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k-E standard 




(d) Radial profile at x/r exit = 51.96 





(c) Radial profile at x/r ex it = 26.93 (f) Radial profile at x/r ex i t = 121.30 

Figure 4. Comparison of computed and measured centerline/radial velocity profiles for Eggers nozzle. 
Computations used SA model, k-e linear model and k-e linear model enhanced with pressure dilatation 
(PD) term and compressibility corrections of Sarkar and Wilcox. 
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(a) Centerline profile 
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(d) Radial profile at x/r ex it = 51.96 






(c) Radial profile at x/r ex it = 26.93 (f) Radial profile at x/r ex i t = 121.30 

Figure 5. Comparison of computed and measured centerline/radial velocity profiles for Eggers nozzle. 
Computations used SA model, SST model and SST model enhanced with pressure dilatation (PD) term 
and compressibility corrections of Sarkar and Wilcox. 
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(a) Mach number 


(b) Stagnation pressure ratio, P t /P t , c 


Figure 6. Comparison of Seiner nozzle’s computed and measured centerline profiles for nozzle plenum 
stagnation temperature of 563 °R. 
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Figure 7. Contours of Seiner nozzle’s computed solution variables on the symmetry plane for nozzle plenum 
stagnation temperature of 2009 °R. 
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(b) Stagnation pressure ratio, P t /Pt, c (b) Stagnation pressure ratio, P t /Pt, c 
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(c) Stagnation temperature ratio, T t /T t c 

Figure 8. Comparison of Seiner nozzle’s computed 
and measured centerline profiles for nozzle 
plenum stagnation temperature of 2009 °R. 
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(c) Stagnation temperature ratio, T t /T u 

Figure 9. Comparison of Seiner nozzle’s computed 
and measured centerline profiles for nozzle 
plenum stagnation temperature of 1359 °R. 
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